Purpose

  • Explore sexual size dimorphism in hummingbird species from Costa Rica

  • Evaluate Rensch’s Rule

 

Report overview

 

Explore data

# read data
dat <- as.data.frame(read_excel("./data/raw/supplementary data.xlsx"))

# read tree
tree <- read.tree("./data/raw/imputed hummer tree 339 spp.tre")

dat$`scientific name`[dat$`scientific name` == "Phaetornis guy"] <- "Phaethornis guy"
dat$`scientific name`[dat$`scientific name` == "Phaetornis longirostris"] <- "Phaethornis longirostris"
dat$`scientific name`[dat$`scientific name` == "Heliomaster longisrostris"] <- "Heliomaster longirostris"
dat$scientific_name <- gsub(" ", "_", dat$`scientific name`)

tree$tip.label[tree$tip.label == "Amazilia_amabilis"] <- "Polyerata_amabilis"
tree$tip.label[tree$tip.label == "Hylocharis_eliciae"] <- "Chlorestes_eliciae"
tree$tip.label[tree$tip.label == "Chlorostilbon_canivetii"] <- "Cynanthus_canivetii"
tree$tip.label[tree$tip.label == "Elvira_cupreiceps"] <- "Microchera_cupreiceps"
tree$tip.label[tree$tip.label == "Calliphlox_bryantae"] <- "Philodice_bryantae"
tree$tip.label[tree$tip.label == "Amazilia_edward"] <- "Saucerottia_edward"
tree$tip.label[tree$tip.label == "Amazilia_saucerottei"] <- "Saucerottia_hoffmanni"
tree$tip.label[tree$tip.label == "Amazilia_saucerottei"] <- "Saucerottia_hoffmanni"
tree$tip.label[tree$tip.label == "Amazilia_candida"] <- "Chlorestes_candida"
tree$tip.label[tree$tip.label == "Elvira_chionura"] <- "Microchera_chonura"

# setdiff(unique(dat$scientific_name), tree$tip.label)

subtree <- drop.tip(phy = tree, tip = setdiff(tree$tip.label, unique(dat$scientific_name)))

names(dat) <- gsub(" ", ".", names(dat))

subtree2 <- subtree
subtree2$tip.label <- gsub("_", " ", subtree2$tip.label)
ggtree(subtree2, layout = "circular", col = viridis(10)[7]) + geom_tiplab(size = 4) +
    theme(plot.margin = unit(c(30, 10, 30, 10), "mm"))

Run statistical models

Size comparison among sexes

log(females) ~ log(males)

(note that in the data supplied ln.SEX.weight variables are in natural log)

# model settings
v.cv.tree <- ape::vcv.phylo(subtree)

prior <- c(prior(normal(0, 10), "b"), prior(normal(0, 50), "Intercept"),
    prior(student_t(3, 0, 20), "sd"), prior(student_t(3, 0, 20), "sigma"))
male_x_mod <- brm(ln.female.weight ~ ln.male.weight + (1 | gr(scientific_name,
    cov = v.cv.tree)), data = dat, family = gaussian(), data2 = list(v.cv.tree = v.cv.tree),
    prior = prior, iter = iter, control = list(adapt_delta = 0.99,
        max_treedepth = 15), file = "./data/processed/male_vs_female_size")
html_summary(read.file = "./data/processed/male_vs_female_size.rds",
    n.posterior = 1000)
male_vs_female_size.rds”
ln.female.weight ~ ln.male.weight + (1 | gr(scientific_name, cov = v.cv.tree))
b_prior sd_prior iterations chains thinning warmup diverg_transitions rhats > 1.05 min_bulk_ESS min_tail_ESS seed
1 normal(0, 10) student_t(3, 0, 20) 10000 4 1 5000 0 0 10742.36 11192.23 821034228
Estimate Rhat Bulk_ESS Tail_ESS CI_low CI_high
Intercept 0.198 1 10742.36 11192.23 0.085 0.316
ln.male.weight 0.838 1 15190.31 13900.16 0.780 0.896

log(males) ~ log(females)

female_x_mod <- brm(ln.male.weight ~ ln.female.weight + (1 | gr(scientific_name,
    cov = v.cv.tree)), data = dat, family = gaussian(), data2 = list(v.cv.tree = v.cv.tree),
    prior = prior, iter = iter, control = list(adapt_delta = 0.99,
        max_treedepth = 15), file = "./data/processed/female_vs_male_size.RDS")
html_summary(read.file = "./data/processed/female_vs_male_size.RDS",
    n.posterior = 1000)
female_vs_male_size.RDS”
ln.male.weight ~ ln.female.weight + (1 | gr(scientific_name, cov = v.cv.tree))
b_prior sd_prior iterations chains thinning warmup diverg_transitions rhats > 1.05 min_bulk_ESS min_tail_ESS seed
1 normal(0, 10) student_t(3, 0, 20) 10000 4 1 5000 0 0 11919.28 13094.22 1392122190
Estimate Rhat Bulk_ESS Tail_ESS CI_low CI_high
Intercept -0.148 1 11919.28 13094.22 -0.291 -0.008
ln.female.weight 1.139 1.001 15511.00 14873.94 1.062 1.218

SSD (Lovich-Gibbons ratio) vs male size

  • SSD represented as aboslute Lovich-Gibbons ratio
  • sex size as log of weight

Absolute SSD vs log male size

male_vs_ssd_mod <- brm(aboluteSSD ~ ln.male.weight + (1 | gr(scientific_name,
    cov = v.cv.tree)), data = dat, family = gaussian(), data2 = list(v.cv.tree = v.cv.tree),
    prior = prior, iter = iter, control = list(adapt_delta = 0.99,
        max_treedepth = 15), s = "./data/processed/male_vs_abs_SSD.RDS")
html_summary(read.file = "./data/processed/male_vs_abs_SSD.RDS", n.posterior = 1000)
male_vs_abs_SSD.RDS”
aboluteSSD ~ ln.male.weight + (1 | gr(scientific_name, cov = v.cv.tree))
b_prior sd_prior iterations chains thinning warmup diverg_transitions rhats > 1.05 min_bulk_ESS min_tail_ESS seed
1 normal(0, 10) student_t(3, 0, 20) 1000 4 1 500 0 0 1486.691 1315.158 631200404
Estimate Rhat Bulk_ESS Tail_ESS CI_low CI_high
Intercept 0.009 1.001 1486.691 1315.158 -0.095 0.111
ln.male.weight 0.065 1 1871.538 1425.567 0.010 0.121

absolute SSD vs log female size

female_vs_ssd_mod <- brm(aboluteSSD ~ ln.female.weight + (1 | gr(scientific_name,
    cov = v.cv.tree)), data = dat, family = gaussian(), data2 = list(v.cv.tree = v.cv.tree),
    prior = prior, iter = iter, control = list(adapt_delta = 0.99,
        max_treedepth = 15), file = "./data/processed/female_vs_abs_SSD.RDS")
html_summary(read.file = "./data/processed/female_vs_abs_SSD.RDS",
    n.posterior = 1000)
female_vs_abs_SSD.RDS”
aboluteSSD ~ ln.female.weight + (1 | gr(scientific_name, cov = v.cv.tree))
b_prior sd_prior iterations chains thinning warmup diverg_transitions rhats > 1.05 min_bulk_ESS min_tail_ESS seed
1 normal(0, 10) student_t(3, 0, 20) 10000 4 1 5000 1 0 18472.36 14149.97 244304415
Estimate Rhat Bulk_ESS Tail_ESS CI_low CI_high
Intercept 0.053 1 18472.36 14149.97 -0.068 0.168
ln.female.weight 0.042 1.001 25317.53 14631.72 -0.027 0.112

 

Takeaways

  • Dimorpism increases with body size in males but not in females

 

 

Session information

## R version 4.1.0 (2021-05-18)
## Platform: x86_64-pc-linux-gnu (64-bit)
## Running under: Ubuntu 20.04.2 LTS
## 
## Matrix products: default
## BLAS:   /usr/lib/x86_64-linux-gnu/atlas/libblas.so.3.10.3
## LAPACK: /usr/lib/x86_64-linux-gnu/atlas/liblapack.so.3.10.3
## 
## locale:
##  [1] LC_CTYPE=pt_BR.UTF-8       LC_NUMERIC=C              
##  [3] LC_TIME=es_CR.UTF-8        LC_COLLATE=pt_BR.UTF-8    
##  [5] LC_MONETARY=es_CR.UTF-8    LC_MESSAGES=pt_BR.UTF-8   
##  [7] LC_PAPER=es_CR.UTF-8       LC_NAME=C                 
##  [9] LC_ADDRESS=C               LC_TELEPHONE=C            
## [11] LC_MEASUREMENT=es_CR.UTF-8 LC_IDENTIFICATION=C       
## 
## attached base packages:
## [1] stats     graphics  grDevices utils     datasets  methods   base     
## 
## other attached packages:
##  [1] viridis_0.6.2       viridisLite_0.4.1   brms_2.18.0        
##  [4] Rcpp_1.0.9          sketchy_1.0.2       brmsish_1.0.0      
##  [7] ggtree_3.2.1        readxl_1.3.1        ape_5.6-2          
## [10] xaringanExtra_0.7.0 rprojroot_2.0.3     formatR_1.11       
## [13] knitr_1.42          kableExtra_1.3.4   
## 
## loaded via a namespace (and not attached):
##   [1] uuid_1.1-0           backports_1.4.1      systemfonts_1.0.4   
##   [4] plyr_1.8.7           igraph_1.3.5         lazyeval_0.2.2      
##   [7] splines_4.1.0        crosstalk_1.2.0      ggplot2_3.4.0       
##  [10] TH.data_1.1-0        rstantools_2.2.0     inline_0.3.19       
##  [13] digest_0.6.31        yulab.utils_0.0.5    htmltools_0.5.4     
##  [16] fansi_1.0.3          magrittr_2.0.3       checkmate_2.1.0     
##  [19] remotes_2.4.2        RcppParallel_5.1.5   matrixStats_0.62.0  
##  [22] xts_0.12.2           sandwich_3.0-1       svglite_2.1.0       
##  [25] prettyunits_1.1.1    colorspace_2.0-3     rvest_1.0.3         
##  [28] ggdist_3.2.0         xfun_0.36            dplyr_1.0.10        
##  [31] callr_3.7.3          crayon_1.5.2         jsonlite_1.8.4      
##  [34] lme4_1.1-27.1        survival_3.2-11      zoo_1.8-11          
##  [37] glue_1.6.2           gtable_0.3.1         emmeans_1.8.1-1     
##  [40] webshot_0.5.4        distributional_0.3.1 pkgbuild_1.4.0      
##  [43] rstan_2.21.7         abind_1.4-5          scales_1.2.1        
##  [46] mvtnorm_1.1-3        DBI_1.1.1            miniUI_0.1.1.1      
##  [49] xtable_1.8-4         gridGraphics_0.5-1   tidytree_0.4.0      
##  [52] stats4_4.1.0         StanHeaders_2.21.0-7 DT_0.26             
##  [55] htmlwidgets_1.5.4    httr_1.4.4           threejs_0.3.3       
##  [58] posterior_1.3.1      ellipsis_0.3.2       pkgconfig_2.0.3     
##  [61] loo_2.4.1.9000       farver_2.1.1         sass_0.4.5          
##  [64] utf8_1.2.2           labeling_0.4.2       ggplotify_0.1.0     
##  [67] tidyselect_1.2.0     rlang_1.0.6          reshape2_1.4.4      
##  [70] later_1.3.0          munsell_0.5.0        cellranger_1.1.0    
##  [73] tools_4.1.0          cachem_1.0.6         cli_3.6.0           
##  [76] generics_0.1.3       ggridges_0.5.4       evaluate_0.20       
##  [79] stringr_1.5.0        fastmap_1.1.0        yaml_2.3.7          
##  [82] processx_3.8.0       purrr_1.0.0          packrat_0.9.0       
##  [85] pbapply_1.6-0        nlme_3.1-152         projpred_2.0.2      
##  [88] mime_0.12            aplot_0.1.7          xml2_1.3.3          
##  [91] compiler_4.1.0       bayesplot_1.9.0      shinythemes_1.2.0   
##  [94] rstudioapi_0.14      gamm4_0.2-6          treeio_1.18.1       
##  [97] tibble_3.1.8         bslib_0.4.2          stringi_1.7.12      
## [100] highr_0.10           ps_1.7.2             Brobdingnag_1.2-9   
## [103] lattice_0.20-44      Matrix_1.5-1         nloptr_1.2.2.2      
## [106] markdown_1.3         shinyjs_2.1.0        tensorA_0.36.2      
## [109] vctrs_0.5.2          pillar_1.8.1         lifecycle_1.0.3     
## [112] jquerylib_0.1.4      bridgesampling_1.1-2 estimability_1.4.1  
## [115] cowplot_1.1.1        httpuv_1.6.6         patchwork_1.1.2     
## [118] R6_2.5.1             promises_1.2.0.1     gridExtra_2.3       
## [121] codetools_0.2-18     boot_1.3-28          colourpicker_1.2.0  
## [124] MASS_7.3-54          gtools_3.9.3         assertthat_0.2.1    
## [127] withr_2.5.0          shinystan_2.6.0      multcomp_1.4-17     
## [130] mgcv_1.8-36          parallel_4.1.0       grid_4.1.0          
## [133] ggfun_0.0.7          minqa_1.2.4          tidyr_1.1.3         
## [136] coda_0.19-4          rmarkdown_2.20       shiny_1.7.3         
## [139] base64enc_0.1-3      dygraphs_1.1.1.6